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ABSTRACT 

Atomic cooling haloes with virial temperatures Tyir > 10^ K are the most plausible sites for 
the formation of the first galaxies and the first intermediate mass black holes. It is therefore im- 
portant to assess whether one can obtain robust results concerning their main properties from 
numerical simulations. A major uncertainty is the presence of turbulence, which is barely re- 
solved in cosmological simulations. We explore the latter both by pursuing high-resolution 
simulations with up to 64 cells per Jeans length and by incorporating a subgrid- scale turbu- 
lence model to account for turbulent pressure and viscosity on unresolved scales. We find 
that the main physical quantities in the halo, in particular the density, temperature and en- 
ergy density profile, are approximately converged. However, the morphologies in the central 
500 AU change significantly with increasing resolution and appear considerably more turbu- 
lent. In a systematic comparison of three diff'erent haloes, we further found that the turbulence 
subgrid- scale model gives rise to more compact central structures, and decreases the amount 
of vorticity. Such compact morphologies may in particular favor the accretion onto the central 
object. 
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1 INTRODUCTION 

The nonlinear interplay of gravity, fluid dynamics, and radiative 
cooling lies at the heart of the continuing challenge to predict the 
properties of the first astrophysical objects of the universe. Scale 
separation, which is one of the most important prerequisites for ro- 
bust analytic or numerical calculations, does not apply: whereas the 
kinetic energy budget is dominated by large-scale modes, fragmen- 
tation and cooling instabilities grow fastest on those scales with the 
strongest density fluctuations, i.e., those close to the Jeans scale. 
In order to capture the relevant degrees of freedom, grid-based hy- 
drodynamical simulations employ adaptive mesh refinement over 
many orders-of-magnitude. 

In addition to the masses and multiplicity of the first stars, 
the question of whether and when supermassive black holes 
(SMBH) can form by direct collapse of a primor dial gas cloud 
belongs to the key topics of current research toh & Haiman' 
2002; Bromm & Loeb 2003; Spaans & Silk 20 06; Begelman et al 
2006; Diikstra et al. 2008; iDiorgovski et all l2008l : IShang et al. 
2010t : .Schleicher et al 20ldnLatif et al.ll201ll) . Massive primor- 



dial haloes assembled at redshifts 10-15 are the plausible sites 
to host the first galaxies formed at the end of cosmic dark ages. 
They are also potential candidates for the formati on of inter- 
mediate mass black holes through direct colla p se (iReesI Il984l: 
Bronim & Loebl l2003l: i Koushia pms et all l2004l: iBegelman et al.l 
20061 : IBegelman & Shlo sman 20091: IVolont"erill20ld) . Thus, their 
understanding is a matter of great astrophysical interest. 



The first stars, on the other hand, are expected to form i n so- 
called minihalos wi th 10^ - 10^ Mo (Abel et al. 2002; Bromm 20041 : 
lYoshida et al.l2008h . In a recent studv. lTurk et al.l(l2012D found that 
an increased resolution per Jeans length results in non-convergence 
of global properties. They discovered that enhanced Jeans reso- 
lution produces higher infall velocities, increased temperatures as 
well as decreased content of molecular hydrogen. As convergence 
has not been achieved even at the highest resolutions, major uncer- 
tainties are present concerning the expected accretion rates, tem- 
perature distributions and the fragmentation behavior. It is there- 
fore important to assess whether similar restrictions apply to other 
systems, in particular the so-called atomic cooling haloes. In the 
study presented here, we therefore assess the convergence behavior 
employing numerical simulations with a resolution of 16, 32 and 
64 cells per Jeans length. We focus on halos exposed to strong Ly- 
man Werner radiation, as these are the candidates for black hole 
formation via direct collapse. We note that these are the high- 
est resolution studies to date, with previous stu dies typically em- 
ploying a re solution of 1 6 cells per Jeans le ngth dWise et al.lf2008l : 
iRegan&Ha ehnelt 20091: IShang et al.ll201Qh . 



The need for high-resolution investigations has been previ- 
ously derived in diff'erent contexts. In particular, Fede rrath et al.l 
(l201lh reported that the turbulent energy as in collapsing gas clouds 
converges only for a resolution of at least 32 cells per Jeans length. 
The turbulent energy in these clouds is released from the gravita- 
tional potential, due to the initial deviations from spherical symme- 
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try (lHovlel[l953l: IScalo & PumD hrevlll982l : iKlessen & Hennebeli3 
l201Ql : lElmegreen & Burk ert 2010). 

It is well known from the field of contemporary star forma- 
tion that turbulence has important implications for the density PDF, 
clump statistics, and angular n iomentum transport in gravitation- 
ally unstabl e ^as clouds (lLarsonlll98ll:lMac Low & KlessenI 



I2OO4I : iMcKee & Ostrikeji2007l : iFederrath & Klessenll2012h ). Dur- 
ing high-redshift structure formation, turbulence was shown to play 
a major role in so-called minihalos. In particular, it regulates the 
angular momentum transport and delays the format ion of a disk 
(lAbel etal.ll2002l :l Brom ml l2004l : lYoihida et al. I [2008b, but also in- 
flue nces the fragmentation behavior Jciark et al.ll201ll : [Smith et aP 
I2OII: Greif et al. 2012). In more massive h alos, the presence of tur - 
bulence was reported bv lGreif et aD (l2008h and lWise et aP (l2008h . 
While the implications were not explored in detail, it is expected 
that it can influence the formation of intermediate mass black 
holes through its impact on disk formation and the angular mo- 
mentum distribution. We note that the pr esence of disks is a cen- 
tral assumption in direct collapse niodels jKoushiappas et al.ll2004l : 



iBegelman et 'al]l2006l : iLodatolzOOVl : IBegelman & Shlosmanll200^ . 
and an efficient means of angular momentum transport is gener- 
ally required to allow t he for mation of a massive central object. 
IBegelman & Shlosmanl (l2009l) therefore invoked the presence of 
nested bar-like instabilities to provide sufficient means for angular 
momentum transport. While our simulations do not exactly support 
this generation mechanism, they show that turbulence is efficiently 
produced during gravitational collapse. 

Such turbulence will not only influence angular momen- 
tum transport, but also amplify existing weak magnetic fields 
via t h e small scale dynamo jKazantseyl 19681: ISubrarnanianl 



1998': Brandenburg & Subramanian' '2005'; Schober et al ' '20121; 
Schleicher et al. 2010; Sur etal. 2010; Federrath et al. 20nj; 



Turk et al boiilSchober et al.l2012h . The enhanced magnetic field 
strength could exert an additional pressure, and further contribute 
to the angular momentum transport. More importantly, this initial 
growth via the small-scale dynamo provides a strong initial field 
on which the ff" - Q dynamo can act to produce the cohe rent fields 
observed today (iBeck et al.lll996l : lArshakian et al.ll2009h . 

We note that an accurate modeling of turbulence not only 
requires high numerical resolution, but also a consistent treat- 
ment of the unresolved scales. For this purpose, we explore the 
implications of a turbul ence subgrid- scale (SOS) model for the 
turbulent kinetic energy ([Maier et al ] |2009l : llapichino etaDl201ll : 
ISchmidt & Federrathll2oT Ih. Its main features are non-ideal terms 
in the fluid equations produced by an efl'ective turbulent viscos- 
ity, an effective turbulent pressure, and fully time and space depen- 
dent dissipation of turbulent kinetic energy into thermal energy. All 
of these terms make order unity contributions for transonic flows, 
which is the case here. In this study, we therefore employ high- 
resolution studies of massive halos to explore their central proper- 
ties. A central question is in particular for which properties con- 
vergence can be achieved, which is directly relevant for the direct 
collapse scenario concerning the formation of intermediate mass 
black holes. 

Our paper is organized as follows. In the next section, we de- 
scribe the simulation setup and the numerical methods employed. 
In the 3rd section of the paper, we present the results obtained in 
this study. In the last section of this article, we summarize our main 
results and confer our conclusions. 



2 NUMERICAL METHODS AND SIMULATION 
DETAILS 

A modified version of the Enzo code including the subgrid- scale 
model for unresolved turbulent fluctuations (see below for a de- 
scription) has been used to perform the simulations presented in 
this work. Enzo is an adaptive mesh refineme nt (AMR), parallel , 
grid-based cosmolo gical hydrodynamics code dO'Shea et al ]l2004l : 
iNorman et al.ll2007l) . It can run on massively parallel systems and 
has been used for a wide variety of astrophysical applications. The 
message passing interface (MPI) is used to achieve portability and 
scalability on different systems. The computational domain is dis- 
cretized into nested grid cells. It has two exchangeable grids, a uni- 
form grid and a block-structured adaptive grid. We use a split hy- 
dro solver with a 3rd order piece- wise parabolic (PPM) method for 
hydrodynamical calculations. The dark matter N-body dynamics 
is solved using the particle-mesh technique. A multigrid Poisson 
solver is employed for the self-gravity computations. 



2.1 Initial conditions 

The simulations are started with the cosmological initial condi- 
tions generated from Gaussian random fields. We employ the inits 
package available with the public version of the Enzo code to cre- 
ate nested grid initial conditions. Our simulations start at redshift 
z = 99 with a top grid resolution of 128^ cells and we select 
the ma ssive halo at redshift 15 using the halo finder of Tur k et all 
(I20T if). Two initial nested levels of refinement are subsequently 
added each with a resolution of 128^ cells. Our simulation box has a 
cosmological size of 1 Mpc h~^ and is centered on the massive halo. 
In total, we initialize 6291456 particles to compute the evolution of 
the dark matter dynamics and have a final dark matter resolution of 
300 Mq. While our dark matter halo is thus well-resolved, we note 
that additional fluctuations could be present in case of a higher res- 
olution in dark matter. The parameters for creating the initial condi- 
tions and the distribution of baryonic and dark matter co mp onent s 
are taken from the WMAP seven years data CJarosik et al.y201ll) . 
We further allow additional 27 levels of refinement in the central 
62 kpc region of the halo during the course of simulation. It gives 
us a total effective resolution of 3 AU in comoving units. The res- 
olution criteria used in these simulations are based on the Jeans 
length, the gas over-density and the particle mass resolution. The 
grid cells matching these requirement are marked for a refinement. 

The simulations conducted in this work mandated the Jeans 
length resolution of 16, 32 and 64 cells throughout their evolu- 
tion. This criterion was applied during the course of simulations 
to ensure that all physical proc esses like shock waves and Truelove 
criterion jTruelove et al.ll 1997b are well resolved. We stop the sim- 
ulations after they reach the maximum refinement level and start to 
violate the Jeans criterion. The results at later stages would not be 
reliable. 



2.2 Chemistry 

To include the primordial non-equilibrium chem- 
istry, the rate equations of the following 9 species: 
H, H+, He, He+, He++, e", H", H2, HJ are self-consistently 
solved in the cosmological simulations. We make use of H2 
photo-dissociating background UV flux implemented in the Enzo 
code. An external UV field of constant strength 10^ in units of J21 
is used in the simulations. We presu me that such flux is g enerated 
from a nearby star forming halo JPiikstra et al.l I2OO8I) and is 
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emitted by Pop III stars with a thermal spectrum of 10^ K. We 
include several cooling and heating mechanisms like collisional 
ionization cooling, radiative recombination cooling, collisional 
excitation cooling, H2 cooling as well as H2 formation heating. 
The chemistry so l ver used in th i s wor k is a modified version of 
lAbel et all ll997h : lAnninos et alJ ( Il997h . 

We note that for columns above 10^^ cm~^, the gas becomes 
optically thick to Lyman a photons, providing a potential compli- 
cation for the further evolution. In fact, ^Spaans & Silk ( 2006) sug- 
gested that the latter may eff'ectively stop the cooling and provide a 
transition to an approximately adiabatic regime. The latter was ex- 
plored in detail by Schleicher et al. ( 2010), finding that at this point, 
additional processes become relevant, including the two-photon de- 
cay (2s- Is transition) and H~ formation cooling. Eff'ectively, the 
temperature evolution is then very close to the evolution obtained 
from optically thin Lyman a cooling. We also note that for a stellar 
spectrum of 10^ K H2 is mainly dissociated by the Solomon process 
while for the stellar spec t rum of 10^ K the main diss ociation route 
is H" dShang et al.ll201Ql : IWolcott-Green et al.ll201lh . In principle, 
one would of course expect the presence of both contributions. The 
details of these processes however would not matter as long as H2 
is efficiently dissociated. 



2.3 SGS turbulence Model 

Due to the high Reynolds numbers relevant for astrophysical sys- 
tems, it is not possible to resolve all scales down to the dissipa- 
tive scale even with adaptive mesh refinement techniques. Turbu- 
lence cascades from coarser grids corresponding to large scales 
down to the center of the structure forming haloes without be- 
ing properly accounted for. In engineering and many other disci- 
plines of computational fluid dynamics subgrid scale turbulence 
models are used to represent the effect of unresolved turbulence 
on resolved scales. The unresolved turbulence has gained lot of in- 
terest in astrophysical simu l ations JScannapieco & Bruggenll2008l : 
lOppenheimer & Davel2009l : lMaier et al.l2009l) . To compute the un- 
resolved turbulence on grid scales in our simulation, we use th e 
subgrid scale (SGS) turbulence model by ISchmidt et all (l2006h . 
This SGS model is based on a mathematically rigorous approach 
separating the resolved and unresolved scales, and connecting them 
via an eddy-viscosity closure for the non-linear energy transfer 
across the grid scale. The turbulent viscosity is given by the grid 
scale and the SGS turbulence energy, i. e., the kinetic energy asso- 
ciated with numerically unresolved turbulent velocity fluctuations. 
The equations for compressible fluid dynamics are decomposed 
into resolved (large scales) and unresolved (small scales) parts us- 
ing the filter formalism proposed by Germano (1992) in terms of 
density weighted quantities. Applying the filtering mechanism to 
the fluid equati ons and solving them in comoving coo rdinate sys- 
tem, we obtain 
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here a is the scale factor and the quantities with ~ are the comov- 
ing gas density, velocity and pressure etc. atj is the viscous stress 
tensor, q is turbulent velocity, X describes the eff'ect of unresolved 
pressure fluctuations, e accounts for the dissipation of kinetic en- 
ergy and 77 is the dynamic viscosity. Equations (1-3) are the mass, 
the momentum and the energy conservation equations. Equation (4) 
solves the evolution of SGS turbulent energy (ct). The quantities D, 
X, 6, r and f (v/, Vj) are unknowns and are computed in terms of 
closure relations, i.e., functions of the filtered flow quantities and 
the turbulent energy. The expressions for the closure terms are the 
following: 
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Further details and numerica l validations of the app lied closures are 

! ;iven in dedicated studies of ISchmidt et al.l (l2006h and lMaier et al.l 
200^ and are beyond the scope of this article. 

For our simulations, the coefficients of the model are cali- 
brated against subsonic compressible turbule nce simulations, com - 
parable to the regime in our simulations fsee lSchmidt et al.|[2006h . 
To apply the SGS model in cosmological AMR simulations, the 
method o f adaptively refined large eddy simulations is employed 
dMaier et al .1120091) . Since SGS turbulence energy depends on the 
grid scale, which varies on adaptive meshes, energy must be ex- 
changed with the numerically resolved velocity field when the grid 
is refined or de-refined. This is achieved by assuming the Kol- 
mogorov two-thirds law for the scaling of the turbulent veloc- 
ity fluctuations. As long as the turbulence is subsonic and nearly 
isotropic locall y, this is a reasonable assumption. In contrast, the 
method used bv lGrav & Scannapiecol (l201lh does not calculate the 
turbulence energy on the grid scale, but the energy associated with 
a characteristic length scale of buoyant bubbles. Our SGS model 
completely neglects gravity on unresolved length scales and as- 
sumes the turbulent cascade from larger resolved scales as the dom- 
inant source of unresolved turbulence (for the eff'ects of buoyancy 
on subgrid scales, see also Schmidt et al. 2006). Of course, this re- 
quires that turbulence is sufficiently resolved. The different resolu- 
tions considered in this study allow us to estimate biases from scale 
separation between the production of turbulence by gravity and the 
grid scale. 
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Figure 1. This figure shows the radially binned spherically averaged radial profiles for the halo C with and without SGS turbulence for diff^erent Jeans 
resolutions. The dashed lines show normal runs while the solid lines represent the cases with SGS turbulence as depicted in the legend. The upper left panel 
of the figure shows the density radial profiles. The temperature radial profiles are depicted the upper right panel. The bottom left panel of the figure shows the 
averaged total energy radial profiles. The averaged radial profiles of SGS energy are shown in the bottom right panel of the figure. 



3 RESULTS 

3.1 Resolution comparison 

In this study, we have performed simulations resolving the Jeans 
length by 16, 32 and 64 cells (hereafter called Ji6, J32. J64 respec- 
tively) and compared the results with SGS turbulence model. The 
properties of the halo for different Jeans resolutions with and with- 
out the SGS turbulence model are shown in figure[T] The results are 
examined for the same peak density (i.e., 10~^^ gcm~^). The aver- 
aged density radial profiles for all the cases are depicted in the top 
left panel of figure[T] The density profiles show almost behavior 
according to the expectation of a isothermal collapse. The bumps 
in the density profile for J16 and J32 cases indicate the presence of 
under-resolved density clumps in the vicinity of the central object. 
The overall density profiles agree well for diflferent resolutions. 

The top right panel of figure [T] shows the averaged tempera- 
ture radial profiles. It can be noticed that for all the cases the halo 
is heated up to its virial temperature and then cools by atomic line 
cooling. The temperature in the center of the halo remains about 
7000 K. This shows that thermal properties of the halo are inde- 
pendent of resolution as well as SGS turbulence. The total energy 
for diflferent Jeans resolutions and SGS turbulence cases is shown 
in the bottom left panel of figure [T] The value of total energy in- 
creases at larger radii and becomes almost constant with a value 



of 3 X 10^^ erg/g in the center of a halo. It is worth noting that 
these quantities are approximately converged, and do not depend 
strongly on the subgrid-scale model. Unlike in minihaloes, the re- 
sults of such systems are thus likely more robust, as the thermal 
evolution is close to isothermal and less sensitive to minor changes 
in the dynamics. 

While the total energy is roughly the same for the three runs 
without SGS model, we see a non-monotonous behavior with res- 
olution if the SGS model is applied. However, by considering the 
morphology of the halo in figure |2l it becomes clear that there are 
only little or no turbulent structure in the lower-resolution runs. For 
J16, turbulence is definitely under-resolved, while J32 could be an 
intermediate case. Consequently, the turbulent cascade cannot be 
sufficiently resolved in these runs. For the highest resolution case 
(i.e., J64), the structure inside the halo appears to be well resolved. 
The central region of a few hundred AU in size has become highly 
turbulent and clumpy. The turn- around in the trend for the total en- 
ergy profiles indicates that Kolmogorov-like turbulence begins to 
be resolved for J64. This is plausible because the scale separation 
between the production of turbulence by gravity and grid scale is 
almost two decades, but the lower decade is strongly aflfected by 
numerical viscosity. 

Our resul ts are also consistent with earlier studies 
(iFederrath et all l201lh . The profiles of the SGS turbulence 



High resolution studies of massive primordial haloes 5 




Figure 2. The figure shows the density projections for the halo C with and without SGS turbulence. The state of simulations for the central 500 AU region of 
the halo is illustrated for various Jeans resolutions. The panels from top to bottom represent the Jeans resolution of 16, 32 and 64 cells respectively. The left 
panels present the normal runs while the cases with SGS turbulence are shown in the right panels. 
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Figure 3. The figure shows the radially binned spherically averaged radial profiles for halos A, B and C with and without SGS turbulence. The solid lines show 
the normal runs while the dashed lines represent the cases with SGS turbulence as shown in the legend. The upper left panel of the figure shows the density 
radial profiles. The temperature radial profiles are depicted in the upper right panel. The bottom left panel of the figure shows the total energy radial profiles. 
The averaged SGS energy radial profiles are depicted in the bottom right panel of the figure. 



energy for the different resolutions are shown in the bottom right 
panel of figure [T] Since the grid scale decreases as the refinement 
levels become higher, the expectation based on the Kolmogorov 
scaling law would be that the SGS turbulence energy decreases 
toward the center and with the maximal resolution. However, 
this applies only to homogeneous turbulence. Since turbulence 
production tends to be stronger in the center of the halo, the 
resolution-dependence is compensated and yields a nearly constant 
SGS energy in the central region. The comparable plateaus of the 
SGS energy for J32 and J64 also indicate that the enhancement of 
SGS turbulence production in the J^a case compensates the smaller 
grid scale in comparison to J32. The SGS energy has a peak around 
radii of 5 x 10^^ cm (2 kpc), which coincides with the drop in the 
temperature profiles. The outward decrease of the SGS energy at 
larger radii shows that turbulence is mainly produced as a result of 
the collapse, but not by hydrodynamical instabilities outside of the 
halo. 

Furthermore, comparisons of the morphology of the halo for 
the three different Jeans resolutions reveal significant differences 
between the runs with and without SGS model (see figure [2]). This 
is most obvious from the low-resolution runs that the SGS model is 
limited by the poorly resolved turbulence. However, we will con- 
centrate on the highest-resolution case for quantitative comparisons 
in the next section. 



3.2 Study of different Haloes 

3. 2. 1 Global Dynamics of collapse 

We have performed six cosmological simulations for three different 
haloes (named A, B and C) with a constant strength J21 = 10^ of 
the H2 photo-dissociating radiation field. The masses of the haloes 
and their collapse redshifts are listed in table [T] The results ob- 
tained from cosmological simulations conducted in this work are 
presented in the following subsections. The density fluctuations 
collapse under the gravitational instability as they decouple from 
the Hubble flow. These small fluctuations merge with each other 
to form larger haloes in accordance with the standard paradigm of 
structure formation. In the early phases of the collapse gas falls in 
the dark matter potentials and gets shock-heated during the non- 
linear evolution phase. Gravitational energy of the halo is continu- 
ously transferred to kinetic energy of the gas and dark matter during 
the course of virialization. 

Here we report the study of three haloes of different masses 
but with same Jeans resolution to compute the variation from halo 
to halo. Again, we compare our results with SGS turbulence model. 
The averaged density radial profiles for the haloes A, B and C with 
and without SGS model, centered at the densest cell are shown in 
the upper left panel of figure [3] The maximum density in our simu- 
lations is a few 10~^^ gcm~^. It can be seen from the figure that all 
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Figure 4. The left panel of the figure shows the H2 abundance for three diff'erent haloes. The solid lines show the cases with no SGS turbulence and the dashed 
lines show the cases with SGS turbulence. The corresponding electron fraction is shown in the right panel. 
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Figure 5. The left panel of this figure shows the density projections for the central 500 AU region of the halo. The images are shown for three diff'erent halos 
A, B and C from top to bottom respectively. The left side of left panel shows the normal runs while the right side depicts the cases with SGS turbulence. The 
right panel of the figure shows the density weighted vorticity projections corresponding to the density projections in the left panel. It can be noticed that halos 
with SGS turbulence model have more compact structures. 
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Figure 6. The properties of the clumps are plotted against their masses in this figure. Diamonds show the data points for normal runs while data points for 
SGS turbulence cases are represented by triangles. Black, red and blue colors of the symbols represent the data for three diff'erent halos (i.e., A, B and C). The 
black solid lines are the fits to normal cases and the dot-dashed red lines are the fits to SGS cases. These clumps have sub-solar masses and are gravitationally 
unbound. 
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Table 1. The halo masses and their collapse redshifts are listed in this table. 



Model 


Mass 


Collapse redshift 




Mo 


z 


A 


8.06 X 10^ 


11.9 


B 


4.3 X 10^ 


11.3 


C 


3.2 X 10^ 


14.1 



cases follow almost behavior which would be expected from a 
isothermal collapse. There is a bump in the density profile of halo 
A with SGS turbulence model which is an indication of fragmen- 
tation. In the very central region the density profile is flat which 
corresponds to the local Jeans length in all cas es. These density 
profiles are comparable to the previous studies (lAbel et aLl[2002l : 
IWise et alJl2008l : lTurk et al1l2009Ll2012h . 

The average radial profile of total energy is shown in the bot- 
tom left panel of figure [3] The total energy of the system is a few 
times 10^^ erg g~^ for all the cases. The figure shows the total en- 
ergy for three diff'erent halos with and without SGS turbulence is 
converged. The increase in the total energy radial profile towards 
the smaller radii is due to gas infall in the center of halo. 

The specific subgrid scale energy for the three diff'erent halos 
is shown in the bottom left panel figure [3] At larger radii, the build- 
up of turbulent energy increases sharply because of the turbulence 
cascade and then gets saturated at smaller radii due to the enhanced 
turbulence dissipation rate. This evolution of SGS e nergy is ac- 
cordi ng to the expectations of large eddy simulations (iMaier et al] 
l2009l) . 

3.2.2 Thermodynamics 

The thermal evolution of the gas for diff'erent halos with same Jeans 
resolution is shown in the right panel of figure [3] During the pro- 
cess of virialization, the gas is heated up to its virial temperature 
(i.e., > IC^K) and subsequently cools by Lyman alpha radiation. 
Consequently, gas collapses almost isothermally with temperatures 
around 8000 K. The temperature profiles for all the cases with 
and without SGS turbulence are similar. There is n o significant 
turbul ent heating for SGS turbulence cases as seen in lMaier et al.1 
(l2009l) due to highly efficient atomic line cooling at these temper- 
atures. The dissimilarities in the temperature profiles at larger radii 
appear due to the diff'erence in the halo masses. The ubiquity of 
intense Lyman Werner radiation photo-dissociates the molecular 
hydrogen and H2 cooling remai ns suppressed. O ur result s are in 
agreement with previous stud ies jBromm & Loeb 2003 ; Wi se et all 
I2OO8I : iLatif et alJl201 iL l2012h and according to the expectation of 
theoretical models. 

The H2 abundance is shown in the left panel of figure lU It 
can be noticed that the H2 fraction increases at lower densities due 
to the rise in electron abundance during the non-linear phase of 
the collapse. At intermediate densities, the H2 abundance becomes 
constant as gas cools, recombines and remains neutral with a con- 
stant temperature around 8000 K. The presence of sharp spikes in 
the H2 fraction is due to the shocks occurring at the central densities 
due to collisional dissociation. In general, the H2 fraction is lower 
than the universal value (i.e., 10~^). Therefore, the contribution of 
H2 cooling in the thermal evolution of the haloes studied here is 
negligible. 

The electron abundance corresponding to the H2 fraction is 



depicted in the right panel of figlU It can be seen that the electron 
abundance is correlated with the H2 fraction. At densities above 1 
cm~^, the electron fraction increases because of virialization shocks 
and then continues to decline as gas becomes neutral. Small wig- 
gles in the central electron fraction are triggered by shocks. 

3.2.3 Halo structure 

The state of the simulations with and without the SGS turbulence 
at the collapse redshift is illustrated by the density projections in 
the left panel of figure [S] Significant changes in the morphology of 
haloes are found in the presence of SGS turbulence in all haloes 
(i.e.. A, B and C). It can be noted that haloes are highly turbulent 
and clumpy in both cases. These eff'ects were not seen in the ear- 
lier studies due to the poor Jeans resolution. Our results confirm 
that one needs to resolve the Jeans length with at least 32 cells 
or higher to capture turb ulent velocity fluctuations Jpederrath et al.! 
I2OI1I : lTurketal.ll2012h . Overall, halos in simulations with SGS 
turbulence are more compact and denser than their counterparts. 
The latter is expected due to the presence of an additional viscosity 
term. As demonstrated above, these haloes have a similar thermal 
evolution but the changes in the morphology arise as a consequence 
of unresolved subgrid scale energy computed via the SGS turbu- 
lence model and show substantial variation from halo to halo. The 
SGS energy is about 10 % of total energy budget. Its effect is par- 
ticularly enhanced on small scales, yielding rather diff'erent mor- 
phologies in the presence of the subgrid-scale model. The structure 
of the halo "A" without SGS turbulence clearly shows that dense 
clumps are very well separated from each other and may lead to 
the formation of a binary in this case. The further evolution of 
the simulations becomes computationally very demanding as the 
Jeans mass keeps decreasing. Here, we stopped our simulations af- 
ter reaching the maximum refinement level. We plan to explore the 
further evolution of at least one halo in a companion paper. 

To determine the presence of turbulent velocity fluctuations, 
we have computed the fluid vorticity (i.e., V x v). Density weighted 
projections of the vorticity squared centered at the densest point are 
depicted in the right panel of figure [5] In the center of the haloes, 
large regions with high values of vorticity indicate the ubiquity of 
high turbulent energy. It is also noted that high values of the vor- 
ticity are correlated with the dense regions of the haloes. These 
vorticity plots further suggest the absence of coherent structures. 
The amount of vorticity is higher compared to the SGS turbulence 
cases because of higher turbulent dissipation rates in SGS turbu- 
lence cases. 



3.2.4 Properties of clumps 

In order to quantify the properties of the clumps fou nd during visual 
inspec tion, we have employed the clump finder of I Williams et al.l 
The properties of the clumps with and without SGS turbu- 
lence model for three different haloes (A, B and C) with power law 
fits are shown in figure (6] The top panel of figure [6] depicts that 
clumps are not gravitationally bound as their masses are smaller 
than the Jeans mass and the number of clumps is generally higher 
in no SGS cases. The ratio of clump mass to Jeans mass shows a 
power law behavior (i.e., M/Mj oc M^^). It is interesting that sim- 
ilar trends have been found in different st udies exploring clumps 
in molecular clouds (iBaneriee et al.ll200"9b . We find that in simu- 
lations with the SGS turbulence model, the clumps have slightly 
higher masses as shown by the fit (red line in figure [6]) and the 
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number of low mass clumps is reduced compared to the no SGS 
cases. This is likely an effect of the turbulent viscosity, which pro- 
vides an additional diffusion mechanism that counteracts the for- 
mation of low-mass clumps. The thermal properties of the clumps 
are depicted in figure [6j showing that the clumps in simulations 
with SGS turbulence have almost the same temperature as those in 
the standard setup. The density of the clumps plotted against mass 
shows a bimodal distribution with clumps sitting at higher densi- 
ties following a linear relation with M. The notable difference is 
that clumps in SGS turbulence case have higher densities and the 
power law sharply drops for lower densities. The latter suggests 
that only the high-mass clumps manage to form in the presence of 
an additional turbulent viscosity, providing a distribution with more 
massive clumps on average. 

The velocity dispersion in the clumps increases with the mass 
and follows a M^^^ power law. Again, it is seen that clumps 
with SGS turbulence have lower values of dispersion velocity but 
roughly follow the same trend. The radii of the clumps comply a 
M^^ growth and clumps with SGS turbulence model have larger 
radii in comparison with no SGS cases. In the last panel of figure[6j 
we show that these clumps are well resolved at least by 10^ cells. 



4 DISCUSSION AND CONCLUSIONS 

We have conducted high resolution cosmological simulations using 
the AMR code Enzo for three different halos with Tvir > lO^K ir- 
radiated by a constant strength of photo-dissociating background 
UV flux. In one set of simulations, \ ye used the subgrid scale 
(SGS) turbulence model pr oposed by I Schmidt et al.l (l2006h and 
ISchmidt & FederrathI JlOllb to compute the kinetic energy of nu- 
merically unresolved turbulence and the associated stresses on re- 
solved length scales. For comparison, we run these simulations also 
without the SGS model. Since a high dynamical range is crucial, 
we use two initial nested grids and insert up to 27 additional re- 
finement levels during the course of simulations, corresponding to 
an effective resolution down to sub AU scales. To investigate res- 
olution effects, we applied refinement at 16, 32 and 64 cells per 
Jeans length. The results from three haloes with different masses 
were examined to study the variation from halo to halo for a fixed 
resolution. The main conclusions from this study are the following: 

• The global properties of the halo, in particular the radial pro- 
files, are converged and can be used as a robust input for direct 
collapse models. 

• Turbulent structures are observed for a Jeans resolution of at 
least > 32 cells. 

• The morphology of the halo and its clump properties are 
strongly influenced by taking into account SGS turbulence and typ- 
ically more compact. 

• The clump properties (i.e., M/Mj, velocity dispersion) show a 
power law behavior against clump masses. 

The gas in the atomic cooling halos is heated up to its virial 
temperature where Lyman alpha cooling comes into play and cools 
the gas down to 8000 K. The presence of an intense Lyman Werner 
UV radiation field of 10^ in units of J21 photo-dissociates the H2 
molecules via the Solomon process. We employed grid resolutions 
of 16, 32 and 64 cells per Jeans length to explore the convergence 
of global properties as well as the local morphology. It is impor- 
tant to note that the radial profiles of density, temperature and total 
energy are approximately converged, implying more robust results 



than previously reported for minihalos jTurk et al.ll2012h . We at- 
tribute the latter to the thermodynamics of these halos, which are 
considerably more robust in the presence of strong H2 photodisso- 
ciation (Schleicher et al. 2O100. In this case, the temperature evolu- 
tion remains very close to isothermal, making it rather insensitive to 
local changes in the dynamics. With such a fixed thermal pressure, 
the resulting evolution during the collapse is therefore considerably 
more robust. We note that SGS turbulence does not have a strong 
impact on thermodynamical properties, again as a result of efficient 
Lyman a cooling. The typical unresolved fraction of the turbulence 
energy is ^ 10% of the total energy. 

Our results demonstrate with the highest resolution simula- 
tions that atomic cooling halos become highly turbulent. We com- 
puted the evolution of three different halos and found that the ra- 
dially averaged properties are very similar for three different halos, 
but the morphology of the haloes varies considerably. The intense 
vorticity inside the haloes demonstrates the absence of coherent 
structures in fully turbulent regions. However, since the non-linear 
coupling between gravity and turbulence converts gravitational po- 
tential energy into kinetic energy, an inertial range of hydrodynami- 
cal turbulence can only exist on length scales smaller than the Jeans 
length. For this reason, a sufficient range of length scales smaller 
than the Jeans length must be resolved to observe a turbulent cas- 
cade with an inertial sub-range of the Kolmogorov type (for low 
compressibility) in simulations. This pushes numerical simulations 
to their limits because an extremely high dynamical range is re- 
quired. In our simulations with 16, 32 and 64 cells per Jeans length, 
we find significant differences concerning the central morphologies 
and the amount of turbulent structures. The latter confirms that tur- 
bulence is only marginally res olved even at the highest resolution 
(see also iFederrath et al.]|201lh . This is also reflected by the radial 
profiles of the turbulent energy in the simulations with SGS model, 
for which no clear convergence trend is found. To verify conver- 
gence, higher resolution > 64 cells per Jeans length should be per- 
formed in the future. 

Based on the notion of large eddy simulations, one would 
naively expect that the application of an SGS model reduces the 
range of length scales that has to be resolved. However, large eddy 
simulations are applicable only if the turbulent cascade is partially 
resolved, i. e., at least by a decade in scale space. On top of that, 
numerical dissipation typically necessitates an additional decade. 
In simulations of turbulent collapsing halos, this corresponds to 
scale separation between the the Jeans length, at which energy is 
injected by gravity, and the grid scale. Apart from that, lack of lo- 
cal isotropy at low resolutions poses a problem because we use 
a constant coefficient for the production of SGS turbulence energy 
by shear. This problem could be addressed, fo r example, by a local - 
ized SGS model with varying coefficients fsee lSchmidt et al.l2006h . 
However, this method has not been applied yet in AMR simulation. 

Overall, it is important to note that, while a subgrid- scale 
model does not yield convergence at low resolution, it certainly 
does improve the solution once a sufficiently high resolution is 
reached, such that its central assumptions are fulfilled. In fact, it 
is clear that direct numerical simulations resolving the turbulence 
over a sufficient range of scales cannot be pursued in the near fu- 
ture. Obtaining robust astrophysical results therefore requires both 
high numerical resolution, as well as subgrid models to account for 
effects still below the grid scale. On the basis of the current simu- 
lations, we can already draw relevant conclusions about the effects 
of numerically unresolved turbulence for the Jeans resolution of 64 
cells. By comparing the results with and without SGS cases, it was 
found that the morphology of the halos obtained in the simulations 
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with SGS model is significantly diff'erent from their counterparts. 
In general, their central gas distributions are denser and more com- 
pact compared to the none-SGS cases. The latter provides an in- 
dication that larger accretion rates onto the central object can be 
expected due to the additional turbulent viscosity. In one case, the 
density structure of the halo shows a bimodal distribution if no SGS 
model is applied, which might result in the formation of a binary. 
Such structural implications need to be addressed employing nu- 
merical methods like sink particles or a pressure floor. It is thus 
worth exploring whether statistical diff'erences can also be obtained 
for gravitationally bound clumps, and how much their accretion is 
enhanced via turbulent viscosity. 

We noticed that structures become compact in the presence of 
SGS turbulence. This may have important implications for the ac- 
cretion of mass to the central object, potentially favoring the higher 
accretion rates. In our simulations, we observed the formation of 
turbulent structures for resolutions of > 32 cells per Jeans length 
and no accretion disk at this stage of the collapse. We do not follow 
the evolution of these haloes for longer dynamical times and there- 
fore cannot make statements about the final fate of these haloes. 
According to theoretical predictions, it is likely that formation of 
disk may take place. Numerical investigations of the direct collapse 
scenario thus need to employ a sufficiently high numerical resolu- 
tion as well as a turbulence subgrid model to determine realistic 
accretion rates. 



ACKNOWLEDGMENTS 

The simulations described in this work were performed us- 
ing the Enzo code, developed by the Laboratory for Computa- 
tional Astrophysics at the University of California in San Diego 
( |http://lca.ucsd.edu| ). We thank Matt Turk, Robi Banerjee, John 
Wise and Enzo developers for helpful discussions. This work was 
supported from the SFB 963 (project A 12) Astrophysical Flow In- 
stabilities and Turbulence. We also acknowledge the funding sup- 
port from German Science Foundation. DRGS thanks for funding 
from the Deutsche Forschungsgemeinschaft (DFG) in the Schwer- 
punktprogramm SPP 1573 Physics of the Interstellar Medium un- 
der grant SCHL 1964/1-1. The simulation results are a nalyzed us- 
ing th e visualization toolkit for astrophysical data YT (iTurk et aP 

Imil). 



REFERENCES 

Abel T., Anninos R, Zhang Y., Norman M. L., 1997, NewA, 2, 
181 

Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93 
Anninos R, Zhang Y, Abel T., Norman M. L., 1997, NewA, 2, 
209 

Arshakian T. G., Beck R., Krause M., SokoM D., 2009, A&A, 
494, 21 

Banerjee R., Vazquez- Semadeni E., Hennebelle P., Klessen R. S., 

2009, MNRAS, 398, 1082 
Beck R., Brandenburg A., Moss D., Shukurov A., Sokoloff" D., 

1996, ARA&A, 34, 155 
Begelman M. C, Shlosman I., 2009, ApJ, 702, L5 
Begelman M. C, Volonteri M., Rees M. J., 2006, MNRAS, 370, 

289 

Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1 
Bromm V., 2004, PASP, 116, 103 



Bromm V., Loeb A., 2003, ApJ, 596, 34 

Clark R C, Glover S. C. O., Smith R. J., Greif T. H., Klessen 
R. S., Bromm V., 2011, Science, 331, 1040 

Dijkstra M., Haiman Z., Mesinger A., Wyithe J. S. B., 2008, MN- 
RAS, 391, 1961 

Djorgovski S. G., Volonteri M., Springel V, Bromm V, Meylan 
G., 2008, ArXiv e-prints 

Elmegreen B. G., Burkert A., 2010, ApJ, 712, 294 

Federrath C, Klessen R. S., 2012, ArXiv e-prints 

Federrath C, Sur S., Schleicher D. R. G., Banerjee R., Klessen 
R.S., 2011, ApJ, 731,62 

Germano M., 1992, Journal of Fluid Mechanics, 238, 325 

Gray W. J., Scannapieco E., 2011, ApJ, 733, 88 

Greif T. H., Bromm V, Clark R C, Glover S. C. O., Smith R. J., 
Klessen R. S., Yoshida N., Springel V, 2012, MNRAS, 424, 399 

Greif T. H., Johnson J. L., Klessen R. S., Bromm V, 2008, MN- 
RAS, 387, 1021 

HoyleF., 1953, ApJ, 118,513 

lapichino L., Schmidt W, Niemeyer J. C, Merklein J., 201 1, MN- 
RAS, 414, 2297 

Jarosik N., Bennett C. L., Dunkley J., Gold B., Greason M. R., 
Halpern M., Hill R. S., Hinshaw G., 2011, ApJS, 192, 14 

Kazantsev A. P., 1968, Soviet Journal of Experimental and Theo- 
retical Physics, 26, 1031 

Klessen R. S., Hennebelle R, 2010, A&A, 520, A17 

Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS, 354, 
292 

Larson R. B., 1981, MNRAS, 194, 809 

Latif M. A., Schleicher D. R. G., Spaans M., 2012, A&A, 540, 
AlOl 

Latif M. A., Schleicher D. R. G., Spaans M., Zaroubi S., 2011, 

A&A, 532, A66 
Latif M. A., Zaroubi S., Spaans M., 2011, MNRAS, 411, 1659 
Lodato G., 2007, Nuovo Cimento Rivista Serie, 30, 293 
Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern 

Physics, 76, 125 

Maier A., lapichino L., Schmidt W, Niemeyer J. C, 2009, ApJ, 
707, 40 

McKee C. F, Ostriker E. C, 2007, ARA&A, 45, 565 

Norman M. L., Bryan G. L., Harkness R., Bordner J., Reynolds 

D., O'Shea B., Wagner R., 2007, ArXiv e-prints 0705.1556 
Oh S. R, Haiman Z., 2002, ApJ, 569, 558 
Oppenheimer B. D., Dave R., 2009, MNRAS, 395, 1875 
O'Shea B. W, Bryan G., Bordner J., Norman M. L., Abel T., 

Harkness R., Kritsuk A., 2004, ArXiv Astrophysics e-prints 

0403044 

Rees M. J., 1984, ARA&A, 22, 471 

Regan J. A., Haehnelt M. G., 2009, MNRAS, 393, 858 

Scalo J. M., Pumphrey W. A., 1982, ApJ, 258, L29 

Scannapieco E., Bruggen M., 2008, ApJ, 686, 927 

Schleicher D. R. G., Banerjee R., Sur S., Arshakian T. G., Klessen 

R. S., Beck R., Spaans M., 2010, A&A, 522, A115 
Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJ, 712, 

L69 

Schmidt W, Federrath C, 2011, A&A, 528, A106 
Schmidt W, Niemeyer J. C, Hillebrandt W, 2006, A&A, 450, 
265 

Schmidt W, Niemeyer J. C, Hillebrandt W, Ropke F K., 2006, 
A&A, 450, 283 

Schober J., Schleicher D., Federrath C, Glover S., Klessen R. S., 
Banerjee R., 2012, ApJ, 754, 99 



12 M. A. Latif et al 



Schober J., Schleicher D., Federrath C, Klessen R., Banerjee R., 

2012, Phys. Rev. E, 85, 026303 
Shang C, Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249 
Smith R. J., Glover S. C. O., Clark R C, Greif T., Klessen R. S., 

2011, MNRAS, 414, 3633 
Spaans M., Silk J., 2006, ApJ, 652, 902 
Subramanian K., 1998, MNRAS, 294, 718 
Sur S., Schleicher D. R. G., Banerjee R., Federrath C, Klessen 

R. S., 2010, ApJ, 721,L134 
Truelove J. K., Klein R. L, McKee C. F, Holliman II J. H., Howell 

L. H., Greenough J. A., 1997, ApJ, 489, L179+ 
Turk M. J., Abel T., O'Shea B., 2009, Science, 325, 601 
Turk M. J., Oishi J. S., Abel T., Bryan G. L., 2012, ApJ, 745, 154 
Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., 

Abel T, Norman M. L., 2011, ApJS, 192, 9 
Volonteri M., 2010, A&A Rev., 18, 279 
WilUams J. R, de Geus E. J., Blitz L., 1994, ApJ, 428, 693 
Wise J. H., Turk M. J., Abel T., 2008, ApJ, 682, 745 
Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 

838 

Yoshida N., Omukai K., Hernquist L., 2008, Science, 321, 669 



